Periodic dynamics of fermionic superfluids in the BCS regime 
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We study the non-equilibrium dynamics of a fermionic superfluid in the BCS limit and in the 
presence of a drive leading to a time dependent chemical potential //(£). We choose a periodic 
driving protocol characterized by a frequency uj and compute the fermion density, the wavefunction 
overlap, and the residual energy of the system at the end of n periods of the drive. We demonstrate 
that the BCS self-consistency condition is crucial in shaping the long-time behavior of the fermions 
subjected to the drive and provide an analytical understanding of the behavior of the fermion density 
nk F (where Jcf is the Fermi momentum) after a drive period and for large uj. We also show that the 
momentum distribution of the excitations generated due to such a drive bears the signature of the 
pairing symmetry and can be used, for example, to distinguish between s- and d-wave superfluids. 
We propose experiments to test our theory. 
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I. INTRODUCTION: 

Ultracold atoms provide us with a useful test bed for 
studying equilibrium and non-equilibrium properties of 
interacting many-body systems. The initial focus in these 
systems has been largely on bosonic atoms; in particu- 
lar, the realization and the study of properties of Bose- 
Einstein condensates (BECs) has been the prime sub- 
ject of investigation in the first few years of experimen- 
tal studies on such systems- 1 . In contrast, studies of 
fermionic atoms have gained momentum much later— 
The main experimental obstacle in studying many-body 
effects in fermionic atoms has been the realization of suf- 
ficiently low temperature so as to obtain a gas of quan- 
tum degenerate fermions with T <Tp ~ h 2 n^ 3 /(^m), 
where m is the mass of the atoms and no is their den- 
sity, Tp is the Fermi temperature of the gas and ks is 
the Boltzman constant. However, recent experiments 
have made significant progress in this direction and it 
has been possible to observe the crossover from classical 
to quantum behavior in fermionic gases 4 . The forma- 
tion of Fermi superfluids, which is an interesting many- 
body phenomenon in its own right 5 , required lower tem- 
perature and stronger interactions. It was soon realized 
that the latter can be achieved by utilizing the Feshbach 
resonance phenomenon which allows for tuning of both 
the strength and the sign of the interaction between the 
fermions. A major hindrance in realizing such strong in- 
teractions for bosonic atoms has been three-body losses; 
in contrast, such losses are minimal for fermionic atoms 
due to the Pauli exclusion principle. This allows for the 
possibility of stable fermionic condensates with strong 
inter-particle interaction which acts as a test bed for 
studying Fermi superfluids and, in particular, the BCS- 
BEC crossover in these systems. Several recent experi- 
ments have verified this phenomenon by numerous mea- 
surements in both the BCS and the BEC side of the 
crossover 6 . 



The dynamical properties of Fermi superfluids have 
also received theoretical and experimental attention in 
the recent past. On the experimental side, there have 
been several studies such as probing the expansion of 
Fermi superfluids after a sudden release of the trap 
potential 7 , measurement of collective excitations of these 
superfluids 8 , measurement of the superfluid gap by radio- 
frequency (RF) spectroscopy 9 , and observation of vor- 
tex dynamics^. On the theoretical side, several studies 
were made to study the equilibrium and near-equilibrium 
properties of these systems. In particular, early stud- 
ies concentrated on understanding the crossover phe- 
nomenon by approaching it from the BCS side 11 . These 
have been later supplemented by inclusion of more so- 
phisticated diagrammatic techniques over the BCS mean- 
field theory 12 , study of the effect of presence of a trap 
potential^ 3 -, inclusion of bosonic molecular degree of 
freedom in the BCS Hamiltonian 14 , and use of quan- 
tum Monte Carlo methods^ 5 -. Later works focussed 
on non-equilibrium aspects of these systems based on 
hydrodynamic approach for studying low-lying collec- 
tive excitations^ 6 , vortex dynamics^ 7 -, quench dynamics 
across a BCS-BEC crossover 18 , and properties of dy- 
namic structure factors of these superfluids in the weak- 
interaction regime 19 . 

Non-equilibrium dynamics of closed quantum systems 
have recently received a lot of theoretical attention due 
to the possibility of realizing such dynamics in ultracold 
atom systems 2 ^ - —. Most of such studies in this direction 
have concentrated on bosonic or spin Hamiltonians real- 
ized by bosonic ultracold atoms in optical lattices^ 9 -"—. 
In particular, experimental realizations of Ising-like spin 
model 34 and Bose-Hubbard model 35 have provided impe- 
tus to such theoretical studies. More recently, concrete 
experiments were carried out on the dynamics of bosons 
near the superfluid-insulator transition 36 . The results of 
such experiments are in qualitative agreement with the- 
oretical studies on such systems 31 . Similar attempts of 
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experimental realization of the Ising model have recently 
been undertaken in trapped ion system o 37 i 38 . However, 
such studies have not been carried out extensively on 
fermionic atoms in the superfluid state. 

In this work, we study the response of a fermionic 
superfluid in the BCS regime to a periodic drive. We 
choose a specific driving protocol which leads to a time- 
dependent periodic chemical potential for the fermions 
characterized by a frequency uj: ja(t) = /io + \i a sin(cjt). 
We note that such periodic drives are known to lead to 
a host of interesting phenomena in quantum systems. 
For example, it has been observed that coherent periodic 
driving in a class of integrable quantum many-body sys- 
tems can give rise to novel quantum phenomena like dy- 
namical many-body freezing, where non-monotonic freez- 
ing behavior (with respect to the driving frequency) is 
observed^. A variant of this phenomenon has also been 
predicted for ultracold bosons in optical lattices 39 . The 
aim of the present work is to study the effect of such a 
drive on superfluid fermions. 

The key results that we obtain from such a study 
are the following. First, we show that the BCS self- 
consistency condition plays a crucial role in shaping the 
response of such superfluids to the periodic drive and 
hence establish that the dynamics of fermionic super- 
fluids will be fundamentally different from those of in- 
tegrable systems such as Ising or Kitaev models which 
can be described by Bogoliubov-like Hamiltonians with- 
out the self-consistency condition. We demonstrate this 
by computing the fermion density (which can be easily 
related to the magnetization of the Ising and Kitaev mod- 
els) which displays oscillatory behavior as a function of 
time for the Ising system and approaches a constant at 
long time for the self-consistent BCS system. We also 
derive an analytical formula for the uj dependence of the 
fermion density rik F (or equivalently magnetization rrik F ) 
at the gap edge (where hp is the Fermi momentum) af- 
ter a complete drive cycle and in the limit of large drive 
frequency. Second, we compute the wavefunction overlap 
(and hence the defect density) and the residual energy of 
the systems after single and multiple cycles of the drive 
and discuss the dependence of these quantities on uj. Fi- 
nally, we compute the momentum distribution of the ex- 
citations created due to the drive at the end of one drive 
cycle and show that such a distribution depends on the 
pairing symmetry of the fermionic superfluid. Thus we 
demonstrate that the dynamic response of these super- 
fluid may prove to be a useful tool for determining its 
pairing symmetry. 

The plan of the rest of paper is as follows. In Sec. 
HH we introduce the model and the corresponding BCS 
mean-field equations and provide explicit expressions for 
the observables that we shall compute. In Sec. IIII1 we 
present numerical results for several observables such as 
the defect density, its momentum distribution, and the 
residual energy at the end of a drive cycle and discuss 
their properties. This is followed by an analytical treat- 
ment of the self-consistent problem in Sec. IIVI where we 



obtain an analytical expression for the uj dependence of 
rrik F at high uj. We provide a discussion of our work and 
suggest possible experiments to test our theory in Sec. IVl 
and provide some calculational details in the Appendix. 

II. FORMALISM 

In this section, we introduce the formalism and define 
the main physical observables which we shall compute 
numerically. The Hamiltonian for a gas of interacting 
ultracold fermions in a shallow square optical lattice, in 
the absence of any drive, is given by 

kcr 
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Here c\^ a represent the annihilation operators for 
fermions of momentum k and spin a = {t?i}- The first 
term represents the kinetic energy of the fermions, and 
the second term the four-fermion interaction energy with 
amplitude g > which represents attractive interaction 
between the fermions. Here ek = — 2J^cos(/^) is the 
band energy spectrum for the fermions, the index i takes 
values x and y for d = 2 or x, y, and z for d = 3, and 
/io is the chemical potential. In the rest of this work, 
we shall assume that the trap potential is slowly- varying 
so that a locally constant chemical potential /io = £k F 
can be used to describe the fermions in the trap. In the 
BCS regime and at zero temperature, the ground state 
of the fermions is a superfluid whose excitations can be 
described by the BdG equations 

F(\c\ ( U ^ \ - ( ( 6k ^ A ( k ) \ ( 

w V ^ ) ~ V A *( k ) - mo) J v ^ J ' 

(2) 

where and are the amplitudes of the particle and 
the hole in a BdG quasiparticle and are related to the 
BCS wavefuntion by 

W) = n( uk+t;kc k ct -k)i )- ( 3 ) 

k 

The pair-potential A(k) depends on the pairing symme- 
try and is given by 

A(k) = Ao, s — wave, 

A(k) = Ao[cos(fc £C ) — cos(fc 2/ )], d x 2_ y 2 — wave. (4) 

In the rest of this work, we shall mostly consider s-wave 
pairing except while discussing momentum distribution 
of the defect density in Sec. \UI\ where we shall discuss 
other pairing symmetries. Our analysis, which will be 
detailed in this section, can be easily generalized to other 
pairing symmetries. For the rest of this work, we set 
h=l. 
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For s-wave pairing, the pair potential satisfies the self- 
consistency relation 

A = g^u&k. (5) 

k 

Eq. [2] and [5] admit the well-known BCS solution 

<>?) = ^(l + (-)^)" 2 . (6) 

We now introduce a time-dependent drive, fi(t) = /x + 
H a sin(o;t), so that /i a , lj <C J. This can be achieved 
in typical experimental systems by introducing an addi- 
tional time-dependent harmonic trap potential which is 
sufficiently broad so as to allow for a uniform fermion 
density. In this regime, the response of the system to the 
drive can be described by the time-dependent Bogoliubov 
de-Gennes equation given by 

B ( u k (t)\ _ /(e k -/*(*)) A(k;t) \ 
"UvkW J ~ { A*(k;t) -(e k -M*)) J 

together with the self-consistency condition which, for 
s-wave pairing, reads 

A(k;t) = A(t) = 9 J2<(tW(t). (8) 

k 

In the rest of this work, we consider the system to 
be in the superfluid ground state at t = t{ with 
(^k(^i), Vk(ti)) = (^ k q , v*^) and study its evolution in the 
presence of the periodic drive till a time tf which corre- 
spond to n cycles of the drive tf — nT = 27rn/cj, where 
n is an integer. 

In order to study such dynamics, we focus on the fol- 
lowing key observables. First, we compute the effective 
magnetization m(t) defined as 

m k (t) = (Mt)\r z \Mt))i 
m(t) = ^m k (t) = ^[l-2Wt)| 2 ], (9) 

k k 

where r z is the Pauli matrix in particle- hole space. The 
observable m(t) shall be equal to m{U) after a full drive 
cycle at T — 2tt/uj both in the impulse (where the wave- 
function does not have time to adjust to the drive) and 
adiabatic limit (where the system remains in the ground 
state of the instantaneous Hamiltonian) . Note that m(t) 
is the magnetization of the Ising or Kitaev models de- 
scribed by BdG-like equations in their fermionic repre- 
sentations sans the self-consistency condition 25,28 . For 
BCS fermions, m(t) can be easily related to the time- 
dependent fermion density rto(t) using the relation 

n (t) = ^n ka = 2^|^ k (t)| 2 = l-m(t). (10) 

kcr k 



Thus m(t) proves to be useful in comparing the behavior 
of integrable Ising and Kitaev models with that of the 
non-integrable self-consistent BCS model. To this end, 
we also define the long time average of m(t) 

1 f nT 

Q = lim — / dtx m(t), (11) 

rwoo nT Jo 

which shall also be used for such comparisons. 

The second quantity which we compute is the wave- 
function overlap. To compute this, we first define the 
amplitudes u^(t) and v^. d (t) which correspond to the 
values of and at time t for adiabatic evolution. 
The ground state of H with \i = ji(t) can be written in 
terms of these quantities as 

\r\t)) = n(< d w+^ d (i)4c f _ k )|o), 

k 

^)[^(*)] = ^(l + <-)fc|^)" 2 ,(12) 

where E(k;t) = ^[ek - fi(t)] 2 + A 2 (t). Note that 
I ^ ad (£/■)) = |^ ad (^)) at the end of any integer number of 
drive cycles. Denoting 

\m) = n («*(*) + MtylcU) io>, (is) 

k 

where u\^(t) and v\^(t) are solutions of Eqs. [7] and we 
define the wavefunction overlap F as 

f = \(r d (t)m))\ 2 

= n f * = n + <* (t) Vk (t)\ 2 .(u) 

k k 

The defect density or the density of excitations generated 
due the dynamics at any instant of time can be written 
in terms of F as 

p d = 5> d (k) (15) 

k 

p d (k) = 1 - F k = K d *(t)v k (t) - vt d *(t)u k (t)\ 2 . 

Note that the defect density identically vanishes for adia- 
batic dynamics and thus provides a suitable measure for 
deviation from adiabaticity. 

Finally, we shall compute the residual energy which is 
the additional energy put in the system due to the drive. 
This is defined as the difference between the energy of the 
system at time t and the adiabatic ground state energy 
and can be written as 

Er(t) = £[fcMt)IM*)l^(*)> 
k 

-«(t)IM*)K(*)>] , 

= i [A 2 (t) + A* 2 (t)-2A 2 ] 

[e k -/*(*)] [m k (t) -<*(*)], (16) 

k 
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FIG. 1: (Color online) Plots of the long time average and 
instantaneous values of the amplitude of the order parameter 
A(t). The left panel shows the plot of the long time averaged 
|A(£)| (averaged over 10 cycles of the drive) with respect to 
the drive frequency uj. The right panel shows variation of the 
instantaneous values of |A(£)| with time t (in units of 2ty/uj) 
for several representative values of uj indicated in the legend. 
In both panels, Ao is indicated by a black dashed horizontal 
line. 

where m| d = 1 — 2|^ d | 2 . Note that the residual energy 
also vanishes for adiabatic dynamics. 

Before closing this section, we note that the BCS self- 
consistency condition imparts dynamics to the order pa- 
rameter A which can be written, using Eqs. [7] and [8] as 

A(t) = ^{A*(t)m(t) + 2^[6 k -/i(t)]<(tK(t)}. 

k 

(17) 

If the time dependence of A is ignored or rendered neg- 
ligible, then the system reverts to an ensemble of de- 
coupled two-level systems in momentum space with con- 
stant gap Ao- In this case, the dynamics is described by 
Landau- Zener-Stiickelberg theory^. As we shall see in 
the next Section, we reach this regime for uj /Ao <C 1; 
however, the behavior of a BCS system differs signifi- 
cantly from that of a bunch of decoupled two-level system 
for moderate uj for which uj/ Aq ~ 1. 



III. NUMERICAL RESULTS 

In this section, we discuss the self-consistent numerical 
evaluation of Eqs. [7] and [8] for d = 2 and subsequent 
computations of m(t), Q, |A(t)|, pd, and E r as defined in 
SecHU We have solved Eq. [71 and [8] for BCS fermions in a 
144 x 144 square optical lattice and having unit hopping 
amplitude (J = 1). The equilibrium gap and chemical 
potential has been taken to be Ao = 0.1 and po = 0.01 
respectively. The periodic drive term has been taken to 
be of the form p a sm(ujt) with p a = 0.1. 

We first consider the plot of the gap amplitude |A| in 
Fig. [TJ The left panel of Fig. [T] shows the average gap 
amplitude as a function of the drive frequency uj after an 
average over 10 cycles. We find that the average value of 
the gap amplitude decreases rapidly with increasing fre- 
quency and keeps fluctuating about |A| ~ 0.2A for large 



FIG. 2: (Color online) Left Panel: Plot of Q as a function of 
uj with averaging carried over 10 drive cycles. Right panel: 
Plot of m(t) as a function of ujt/2it for representative values 
of uj indicated in the legend. A few representative values of 
the adiabatic magnetization m adb (t) is shown using crosses. 
In both panels, the initial values of Q and m is indicated by a 
magenta dashed horizontal line and all parameters are same 
as in Fig. [1] 
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FIG. 3: (Color online) Same as in Fig. [2] but for the non-self- 
consistent dynamics. 



^/Ao > 2. The right panel shows a plot of |A(t)|/Ao as 
a function of ujt/2ir. We find that at small uj <C Ao, 
|A(t)| displays oscillatory behavior with maximum and 
minimal values of Ao and 0.9Ao respectively. However, 
for uj > Ao, the behavior of |A(£)| is qualitatively differ- 
ent; it decreases rapidly to near-zero values within the 
first couple of drive cycles {yjt/2i\ < 2) and continues to 
fluctuate around this value for longer drive times, never 
returning close to its original value A(0). We note that 
such a behavior of | A{t) | clearly reflects the importance of 
the self-consistency condition in the dynamics; any anal- 
ysis with |A(t)| ^ Ao at all times is expected to produce 
qualitatively wrong results for uj > A$. 

Next, we plot the effective magnetization m(t) as a 
function of time t and its time average Q as a function of 
the drive frequency uj. For the self- consistent dynamics 
appropriate for fermions in the BCS regime, as shown in 
Fig. [2j there are clearly three regimes, one crossing over 
to the other as uj is increased. For uj <C Ao, m(t) (right 
panel) oscillates with large amplitude, following the drive 
almost adiabatically, resulting in Q = m(0). The oscilla- 
tions, though large, respects the symmetry of the drive. 
As uj approaches Ao, this symmetry is destroyed, result- 
ing Q m(0). For uj > Ao, the oscillatory behavior of 
m(t) with large amplitude is replaced by relaxation to 
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an approximately constant value (with negligible fluctu- 
ations) within few initial cycles (right panel). This con- 
stant value determines the value of Q (left panel), and we 
find that it deviates steadily from mo as uj is increased 
for uj < 2.5Ao. We note that for uj > Ao, the mixing of 
the k modes of the quasiparticle excitations, which origi- 
nates from the presence of the self-consistency condition 
and is therefore absent in Ising or Kitaev systems, is at 
the heart of such a deviation. The freezing of magnetiza- 
tion breaking the symmetry of the driving, seen here for 
uj < Aq, was observed in Ref. 28 for periodically driven 
transverse Ising chain for large amplitude and frequency. 
For very large uj, we of course see the behavior of m(t) 
crossing over to a regime where Q — >• m(0) again - here 
uj becomes too large for the system to react at all, and 
m(t) remains frozen around m(0) for all time (the regime 
sets in beyond uj > 2.4). 

The above behavior is to be contrasted with the non- 
self- consistent case summarized in Fig. [3l Here m(t) al- 
ways executes a large, almost synchronized oscillation, 
approximately following the adiabatic path (black crosses 
in the right panel of Fig.[3j). Naturally, the resulting val- 
ues of Q are close to ra(0) (albeit with some small fluctu- 
ations). This suggests that the synchronous oscillation is 
simply a manifestation of the near- adiabatic nature of the 
dynamics. Synchronization could also occur due to self 
dephasing of the system, after all the transients (some 
of them having power-law tails) have died down, due 
to quantum interference between the modes 2 -^i. But 
such synchronization would appear only in the ujt —> oo 
limit, unlike in the present case, where the effect is visible 
from the very first cycle. The qualitative departure from 
this behavior in the self-consistent case seems to stem 
form the non-adiabaticity induced by the self-consistency 
condition (Eq. [8]) which makes the effective Hamiltonian 
non-linear. The overlapping eigenfunctions of the non- 
linear Hamiltonian makes the criteria for adiabatic be- 
havior much more restricted compared to a linear case 
(see. e.g., Ref. |42| ). We shall address the behavior of 
rrik F (t) in a more quantitative manner in Sec. [IVl where 
we shall show that the value of rrik F after one drive cycle 
decays as 1/uj for uj ^> Ao. 

Next, we consider the behavior of the defect density pd 
as shown in Fig. 2J Here, as expected, the defect density 
becomes significant only for uj > Ao . The plot of the self- 
consistent defect dynamics shown in the left panel of Fig. 
H] demonstrates that the defect density is an oscillatory 
function of uj. The time- averaged defect density shown in 
the right panel of Fig. [4] for both the non-self-consistent 
and the self-consistent dynamics shows that these quan- 
tities display qualitatively similar behavior. A similar 
behavior is seen for the residual energies as can be seen 
from Fig. [5] We find that E r (t) vanishes for uj <C Ao and 
displays oscillatory behavior for uj > Ao- The behavior 
of residual energy and the defect density clearly shows 
that the system wavefunction never comes back to itself 
for any uj; thus BCS superfluids do not seem to support 
dynamic freezing as predicted for superfluid bosons in 
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FIG. 4: (Color online) Left panel: Plot of the instantaneous 
defect density pd(t) as a function of ujt/{2n). The values of 
the drive frequency uj is indicated in the inset. Right panel: 
Plots of the long time average (averaged over 10 drive cycles) 
of the defect density (both the non-self-consistent and the 
self- consistent cases) as a function of the drive frequency uj. 
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FIG. 5: (Color online) Same as in Fig. [H but for residual 
energy E r . 
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Finally, we consider the momentum distribution of the 
defect density at the end of a drive cycle and compare 
such plots for d-wave and s-wave pairing symmetries. 
The generalization of our calculation for d-wave pair- 
ing symmetry is straightforward and constitutes changing 
A -» A k = Aolcos^) - cos(k y )\ in Eqs. [3 and [3 The 
rest of the computation follows exactly as charted out in 
Sec.HH We expect the momentum distribution of the de- 
fect density to be qualitatively different for s- and d-wave 
pairing symmetries. For s-wave, the defect density has a 
uniform pattern in the Brillouin zone as expected from 
the momentum independence of the order parameter. In 
contrast, for the d-wave pairing symmetry, A(k) vanishes 
around k x = ±k y . It is easy to see from Eq. [7] that for 
such momenta, one has 



UkM = ^(^-|k|)e"'^ [ek -^ dt/ , 
Vk(t) = 0(\k\ -k F )e lS ^-^ t,)]dt \ 



(18) 



where u\^(U) = 6(kF — |k|) and v\^(U) = 0(|k| — kp) 
are obtained by solving the BCS equations (Eq. [2j) for 
A(k) = 0, and kp denotes the instantaneous Fermi 
momentum at t = t{. We note that this also implies 
that at the end of a cycle, at t = tf = uj/2tt where 
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fjb(tf) = the phase integrals vanish and one obtains 

u\i(t f ) = 6(k F - |k|) and v w (t f ) = 0(|k| - Thus the 
wavefunction overlap for A(k) =0 at the end of a drive 
cycle becomes unity leading to vanishing defect density 
at the nodes. However, away from the nodes, where the 
A(k) is finite, we expect high density of quasiparticle 
excitations. This qualitative consideration matches well 
with the numerical results shown in the right panels of 
Fig. [6] In contrast, the s-wave pairing symmetry has a 
constant A(k) = Aq and hence leads to a uniform defect 
density pattern as shown in left panels of Fig. [6] This 
results in a qualitative difference between the defect den- 
sity patterns originating from superfluids with the two 
pairing symmetries. We note that although we have ex- 
plicitly calculated the defect density for s- and d x i_ y i- 
wave symmetries in the present work, the approach can 
be straightforwardly generalized to other pairing symme- 
tries. In general, we expect the defect density to display 
a minimum at the position of the node of the gap. Thus 
the momentum distribution of the defect density of a 
Fermi superfluid bears the signature of the positions of 
the nodes of the order parameters on the Fermi surface 
and hence can be used to distinguish between various 
order parameter symmetries. 



IV. ANALYTICAL COMPUTATION OF THE 
MAGNETIZATION 

In this section, we obtain an analytical understanding 
of the behavior of the magnetization or fermion density 
at the gap edge, i.e. at |k| = kp-, after a drive cycle 
in the high-frequency limit. We note that from Refs. 
Eol and |43|, we can conclude that the non self-consistent 
dynamics for a two- level system described in Sec. [II] is 
affected by two phenomena, Landau- Zener tunneling and 
the Stiickelberg phase. In what follows, we derive an 
analogous picture for the s-wave BCS fermions with the 
self- consistency condition. The calculation is carried out 
here for s-wave superfluids but can be easily generalized 
to other pairing symmetries. 

We begin with Eqs. [3and[8j yielding 



*k(*) 
A(t) 



-i [e k - Uk(t) - iA(t)vk(t), 
i[e\<- v(t)]v k (t) -iA*(tK(t), 



(19) 



Defining the terms 

s k (t) = ex p{"^ dt'[e k -/i(t / )]|, 

Uk(t) = a k (t)s k (t), v k (t) = bk^Sfe 1 ^), 

the dynamics of the system can be rewritten as 

a k = -iA(t)b k (t)s^ 2 (t), 
6 k = -iA*(t)a k (t)si(t), 



(20) 
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FIG. 6: (Color online) Top (Middle) Panels: Plot of the mo- 
mentum distribution of the defect density pd(k) after a full 
drive cycle for s-wave [left panels] and d-wave [right panels] 
pairing symmetries. The drive frequencies are specified in the 
inset of the right column. In these plots, lighter shades rep- 
resent higher quasiparticle excitation densities. The bottom 
panel indicates the momentum distribution of the defect den- 
sity (black dotted line for s-wave and red solid line for d-wave) 
along the Fermi surface in the first quadrant (kx , ky > 0) as 
a function of k x = k^ clearly demonstrating the dip in the 
momentum distribution at the node for d-wave pairing. Anal- 
ogous dips occur at other three quadrants at the position of 
the other nodes of A(k). 



which leads to two decoupled second-order differential 
equations for a^(t) and b\^(t) given by 



a k - J 2i [e k - fi(t)] + \ d k + |A(t)| 2 a k = 0, 



k+{2i [e k - n(t)] 



A*(t) 
A*(t) 



6 k +|A(t)| 2 6 k = 0. 



(22) 



(21) The self-consistency condition can be written in terms of 
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Value at t = 0— 


First derivative at t = — 


Second derivative at t = — 


n k = ^ 


= (ik^k q + A v® q ) 


iik = — [E 2 (k) — iji a ui^ 




i>k = +i (/kC - A < q ) 


ii k = - [£ 2 (k) + i/i a u;] < q 


S k = 1 


Sk = — i/k 


S'k = - (/k - VaW) 


6k = 


6 k = -iA n^ q 


6 k = -A < q (A + 2/ k ) 


A = A 


A = 


A = —2iA fi a uj 



TABLE I: Values of system dynamical variables Uk(t), v^(t), Sk(t), b^(t), and A(£), their time derivatives and second derivatives 
at t = where the avoided crossing takes place for adiabatic initial conditions. See Eqs. [6] [8] and [TJ] for details. 



<2k(£) and &k(£) as 



A(t) =#exp 



^0 J u 



;(t)6 k (t)e 2 ^^, (23) 



where /k = ^k — /io and m'(^) = /iasincjt. The initial 
conditions for ak and 6k can be easily obtained from those 
of ?ik and ^k as discussed in Sec. HU 

To obtain an analytical insight into the solution of 
these equations, we note that there is an avoided cross- 
ing at tik = t\ = arcsin(/k/Ma)/^ and that t\ approaches 
zero for large uj and for all momenta on the Fermi sur- 
face. Further if ujt\ <C 1, a condition which is exactly 
satisfied at |k| = hp, we may use the Zener approxima- 
tion fi'(t) ~ [L a ujt for n(t) close to t = t\ when the system 
traverses the avoided crossing 43 45 , and restrict ourselves 
within the adiabatic impulse model where all excitations 
away from the avoided crossing are ignored 4 ^. From the 
definitions in Eq. [20] and Eq. [T9l we can simplify Eq. [21] 
to yield 



-zA(t)6 k (t)e 2i ( /k£ -^ aU;£2 ). 



(24) 



We now assume that uj ^> Aq and define 



X k ty/fJL a LJ - / k / y/ll a U 

0k(*) = 6k(*)A(t)/A =^(t)5k(t)A(t)/A . (25) 
We can use these definitions to simplify Eq. [24] yielding 



z k .A 6> k (^k) Hi- _ 
— — = —i e^g^e 



(26) 



where denotes the amplitude ak after n passages 

across the avoided crossings with A^ being the initial 
amplitude of a k 4Ql43 . 

The integral in the right side of Eq.[23can be evaluated 
by contour integration whose details are charted out in 
the Appendix. This yields 



A 



(i) 




n=0 



n\ (4/i a o;) n d 2n t 



(28) 



t=t k 



where we have used ak(O) 



eq 



from Eq. EUJ and tk = 



/k/ (/^a^)- Note that, in general, £k 7^ and so evaluating 
the modified Landau Zener probability for an arbitrary 
momentum will require knowledge of the system at times 
tk- However, vanishes exactly on the Fermi surface 
(characterized by /k = 0) and can be set to zero for all 
k that lie within 0(fi a uj /vp) around the Fermi surface, 
where vp is the Fermi velocity. In the rest of this section, 
we shall restrict ourselves to this limit. 

The fermion density = 2|S^| 2 = 2|^k| 2 after one 
passage across the avoided crossing can be obtained in 
terms of the modified Landau Zener probability 



(1)|2 _ 



K q ) 2 + Xo|c k | S 



-2«k q Vxo x Re (ci^ 



(29) 



where we have denned 



Thus, the amplitude a k (t) after the system traverses an 
avoided crossing is approximately given by 



A 



(i) 



A 



(o) 




dxk#k(x k )e-^, (27) 



Ck 



^k 



1 Hf 9 



2nz 



n=0 

JL 



n!(4/i a o;) n d 2n t 



t=tk 



7T 

4' 



Xo 



ttA^ 



(30) 
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FIG. 7: (Color Online) Numerical plots (blue crosses) of mk F 
as a function of uj. The system parameters are the same as 
in Fig. [U The red solid line indicates the analytical result 
obtained in Eq. 1391 



Noting that for uj ^> /k, ^k ~ tt/4, and approximating 

(31) 



~ i (-»)" d 2n e k 



n=0 



n! (4/i a cj) n <9 2n £ 



£=0 



one finally gets an expression for the modified Landau- 
Zener probability 



« q ) 2 



Xo|c k | 2 + < q V2xoRe [(1 + i) c k ] 

(32) 



Each of the terms in the sum of Eq. [31] can be obtained 
from table U and higher order derivatives thereof using 
Eq. H3 and either Eq. ED or Eq. |22] at t = 0. We now 
define Pk to be the traditional Landau Zener probability 
for the non self-consistent case 43 45 viz. 



Now, we can write |B, 



(1)|2 



(33) 



(P k /2)e 7k , where 

7k = Xo + m{2« q ) 2 -2 X o|c k | 2 
+2< q y2^Re[(l + i)c k ]} 



(34) 



for large uj. We now investigate regions close to the Fermi 
surface by simplifying Eq. [3TJ retaining only the terms up 
to n = 1 in the expansion. This yields 



2tt 



1 



2/k 
A 



(35) 



where we have used expressions from table HI Taking the 
approximation for Ck in eq[35j substituting its value into 



Eq. M 

yields 



and retaining lowest contributing orders of xo 



IB 



(1)|2 



CI 2 



V2 < q 



(36) 



The occupation amplitude B k "^ is realized after the first 
passage across the avoided crossing and when the sec- 
ond passage begins. The passage starts when the adia- 
batic energy equals the BCS gap Aq. The adiabatic ener- 



gies are given by E(k;t) = ±y [/ k - fi a sin ujt] + |A | 2 , 
which is E(k) (Eq. [6]) with /i replaced by /io + Ma sincjt. 
The passage ends when the velocity of the adiabatic en- 
ergy vanishes, i.e when E(k;t) = at t = 7r/co> or half a 

period. Thus, |B^| 2 is the fermion density after half a 
drive cycle. After one complete period i.e two passages 
across the avoided crossing, the fermion density in the 
adiabatic impulse limit is given by 40 



IB, 



(2) 12 _ 



^0 

(u e Xq) 2 ^ _ 



sin 2 [$ s 



(h 



,e q|2 



CI 2 




^ , (37) 



where in the second line, we have used the fact that 
the Stiickelberg phase <£ st — » 7r/4 for large uj 40 . Thus 
the magnetization ra k ; after one period (or two passages 
across the avoided crossing) can be evaluated using Eqs. 
[9]and[37J yielding 



m 



(2) 



1-41 



f)( 



(38) 



where m k q , the equilibrium magnetization, is given by 
Thus at |k| = kp-, where all approxi- 



rru 



1 



2|<T. 



mations used to arrive at this result are clearly valid, one 



obtains, using u 



eq eq 

1 1 1 



= 1/V2, 



,( 2 ) 



Wife* 



ra - 



m 



ttAq 

2/i a ' 



(39) 



We note rrik F does not depend on the position of k on 
the Fermi surface. This is a consequence of the s-wave 
symmetry of the superfluid order parameter and is not 
going to be present for other pairing symmetries. 

Thus, we find that the frequency dependence of the 
magnetization (or equivalently the fermion density) is 
rrik F ~ uj~ 1 . The magnetization drops off at the same 
manner as the non self-consistent case {i.e the driven 
Ising model where m ~ uj -1 can be obtained from the 
Landau Zener formula). The constant mo which decides 
the rate of the decrease of rrik F with cj, however, is differ- 
ent in the two cases. In the non self-consistent case, A is 
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always constant and so c k = 1/ y/2 exactly as yielded by 
Eq. [31] Thus, « 2%o to lowest order which is four 

times its value for the self-consistent case. We note here 
that although we have concentrated on m^ F , our results 
are expected to be accurate for all m k for which t\ ~ 
and /k <C /i a . The generalization of this treatment to 
n periods seems to be difficult due to the necessity of 
taking into account multiple Stuckelberg phases and we 
leave this issue for a possible future study. 

To check the accuracy of the analytical result, we com- 
pare it (Eq. [39]) with numerical results for the magneti- 
zation on the Fermi surface after one drive period for 
a 144 x 144 square optical lattice which is very close to 
half filling (jao = 0.01) with /i a = Ao = 0.1. At each time 
step of the numerics, the Fermi surface was strobed with 
a tolerance ~ fi a uj. In the limit of /k <C /i a , the variation 
of rrik F is expected to be small within this region of the 
Brillouin zone and an average over all momenta inside 
this region should yield values close to that predicted by 
Eq. [39j The agreement, as shown in Fig. [7J is quite good 
for uj > 3Ao but poor for smaller uj /Ao where some of the 
approximations made in this section are clearly violated. 



V. DISCUSSION 

Experimental verification of our work will require a 
similar momentum distribution measurement as done re- 
cently for fermions on a honeycomb lattice in Ref. |46|. 
A comparison of momentum distribution of the super- 
fluid fermions before and the after the dynamics could 
be used to measure the momentum distribution of the 
defects generated during the drive. Our theory predicts 
that this momentum distribution would depend on the 
pairing symmetry of the superfluid and its pattern would 
be qualitatively similar to that showing in Fig. [6] for s— 
and d x 2_ y 2— wave superfluids. 

In conclusion, we have studied the non-equilibrium pe- 
riodic dynamics of Fermi superfluids in the BCS regime 
within a self-consistent mean-field theory. We have 
shown that proper incorporation of the self-consistency 
condition is crucial for understanding the dynamic prop- 
erties of such systems. This is particulary highlighted by 
studying the behavior of the effective magnetization m(t) 
(or equivalently no) which shows qualitatively different 
behavior for Fermi superfluids (which obey self-consistent 
BCS equations) and Ising or Kitaev spin models (whose 
properties are governed by BCS-like equations without 
the self-consistency condition). We have also studied the 
behavior of defect density, it's momentum distribution, 
and the residual energy for such dynamics. In particu- 
lar, we find that the momentum distribution of the de- 
fect density bears a signature of the pairing symmetry of 
such superfluids. Finally, we have provided an analytical 
derivation of the frequency dependence m^ F at the end 
of one drive cycle in the limit of large drive frequency 
and have shown that rrik F ~ 1/cj for uo ^> Ao- 
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FIG. 8: Contour C = Ci + C2 + C3 + C4 used for evaluating 
the integral in Eq. (|27|) in the complex x-plane. The path 
traversed is indicated by arrows, and the contour consists of 
two isosceles right angled triangles of height xo laid out as 
shown above. 
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Appendix 

Here we provide details of the evaluation of Eq. [271 in 
Sec. [IV[ The integral in the right side of Eq. [23 can be 
evaluated by following the contour C = Ci -f C2 + C3 + C4 
in the complex plane as shown in Fig. [51 Along the 
path Ci, we have dz = dx and exp (—iz 2 ) = exp (— ix 2 ). 
Along the path C2, we have dz = —idy and exp (—iz 2 ) = 
exp [—i(xo - iy) 2 } = exp [—i(x% - y 2 )] x exp (-2xoy), an 
integrand that vanishes when xq — > 00. Along the 
path C3, we have dz = (1 — i)dx and exp(— iz 2 ) = 
exp(— 2x 2 ). Finally, the integrand vanishes along path 
C4 in a way similar to that along C2. Since the contour 
C does not enclose any poles, Cauchy's theorem yields 
§ c dzd\Sz)e~ %z = 0. Thus, taking the limit xq —> 00, 



dx k (9 k (x k )e lx * = 

/oo 
dx k k [(l-i)x k ]e- 2 < (A.l) 
-OO 

Performing a Taylor expansion of # k , 



n=0 



(A.2) 
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and substituting this into the right side of Eq. IA.1I after 
the transformation — >• (1 — i)x\^ each term in the sum 
can be evaluated using Gaussian integrals, yielding 



/oo 
dx k 6 k [(l-i)a; k ]e- 2 ^ = 
-OO 



2 ^ n! (4/i a w) n d 2 ^ 

n=0 v 7 



(A.3) 



t=t k 



Using the above result, we evaluate Eq. IA.1I and hence 
Eq. E3 This yields Eq. [28] used in Sec. Ml 
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